Single-cell combined bioinformatics analysis: construction of immune cluster and risk prognostic model in kidney renal clear cells based on CD8+ T cell-associated genes

Background Kidney cancer is an immunogenic solid tumor, characterized by high tumor burden and infiltration of CD8+ T cells. Although immunotherapy targeting the PD1/CTLA-4 axis has demonstrated excellent clinical efficacy, clinical outcomes in most patients are poor. Methods We used the RNA sequencing data from the GEO database for KIRC GSE121636 and normal kidney tissue GSE131685, and performed single-cell analysis for cluster identification, pathway enrichment, and CD8+ T cell-associated gene identification. Subsequently, the significance of different CD8+ T-cell associated gene subtypes was elucidated by consensus clustering, pathway analysis, mutated gene analysis, and KIRC immune microenvironment analysis in the TCGA–KIRC disease cohort. Single gene analysis identified LAG3 as the most critical CD8+ T-cell-associated gene and its function was verified by cell phenotype and immunohistochemistry in KIRC. Results In the present study, CD8+ T-cell associated genes in KIRC were screened, including GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 and an immunological risk prognostic model was constructed (risk score =  − 0.291858656434841*GZMK − 0.192758342489394*FXYD2 + 0.625023643446193*LAG3 + 0.161324477181591*RGS1 − 0.380169045328895*DUSP4 − 0.107221347575037*TRBV20-1). LAG3 was identified and proved as the most critical CD8+ T cell-associated gene in KIRC. Conclusion We proposed and constructed an immunological risk prognostic model for CD8+ T cell-associated genes and identified LAG3 as a pivotal gene for KIRC progression and CD8+ T-cell infiltration. The model comprehensively explained the immune microenvironment and provided novel immune-related therapeutic targets and biomarkers in KIRC. Supplementary Information The online version contains supplementary material available at 10.1186/s40001-024-01689-8.


Introduction
Kidney cancer is a common malignancy of the urinary system that originates from the renal tubular epithelium.The most common histological subtype of kidney cancer is renal cell carcinoma (RCC), accounting for approximately 90% [1,2], including kidney renal clear cell carcinoma (KIRC, 70%), papillary RCC (10-15%), and chromophobe RCC (5%) [2].The Global Cancer Statistics (2020) stated that RCC accounted for 2.2% of newly diagnosed cancers annually, of which 25 to 30% of patients were diagnosed as advanced or metastatic with a 5-year survival rate of 10%, and 20% to 30% of patients having a propensity for recurrence and metastasis after local operation [3][4][5].
Although the application of immune checkpoint inhibitors (ICIs) as tumor therapeutics has led to major improvements in the RCC systemic treatment, most patients fail to achieve a durable complete response (CR).This could be because RCC is significantly different from other solid tumors in immunogenic features, has a high mutational burden and CD8 + T cell infiltration, and is associated with poor prognosis [10,11].CD8 + T cells constitute the major anti-tumor effective cells in the tumor microenvironment (TME) and exert cytotoxic effects.However, their function is impaired by immunosuppressive cells or molecules in the TME [12].Up-regulation of co-suppressor molecules, including PD1 and CTLA-4, on the surface of CD8 + T cells bind to relevant ligands, ultimately causing CD8 + T cell exhaustion [13,14].Therefore, blocking PD-1-mediated inhibitory signaling by monoclonal antibodies could reverse CD8 + T cell exhaustion, thereby hindering tumor progression.However, this contradicted CD8 + T cell characteristics in RCC [10,11,15].Recent studies have reported that the timing of PD-1 inhibition could negatively affect T-cell priming and memory CD8 + T cell formation, thus contributing to more appropriate timings in RCC immunotherapy [16,17].
No relevant studies are available to explore the function of CD8 + T cell-associated gene sets in the KIRC and its TME.The high accuracy and specificity of single-cell RNA sequencing allows analysis of the functional status of CD8 + T cells and the expression of its associated genes using the single-cell sequencing data of immuneinfiltrated KIRC (GSE121636 [18]) and normal kidney (GSE131685 [19]).In the present study, we developed a CD8 + T cell-associated immunological risk prognostic model using CD8 + T cell-associated marker genes obtained from single-cell sequencing analysis to predict the survival status, tumor immune microenvironment, and immunotherapy responsiveness of KIRC patients, thereby providing a potential target and predictive evidence for immunotherapy.

Single-cell data filtration, pre-processing, and cluster identification
The single-cell sequencing data of KIRC GSE121636 [18] and normal kidney GSE131685 [19] were screened using the Gene Expression Omnibus (GEO) database (https:// www.ncbi.nlm.nih.gov/ geo/).See Table 1 for details.The Seurat package was used to generate objects and filter out poor-quality cells.The standard data pre-processing processes were performed and percentages of gene numbers, cell counts, and mitochondrial sequencing counts were calculated.The filtering criteria were genes with less than only three cells detected and disregarded cells with less than 200 detected gene numbers.Cells with less than 200 or more than 2500 genes detected and those with high mitochondrial content (> 10%) were filtered out as well.We scaled the UMI counts using a scale factor of 10,000 to normalize the library size effect of each cell.Following the log transformation of data, other factors, including "percent.mt," "nCount_RNA, " and "nFeature_RNA" were corrected for variation regression using the Scale-Data function in Seurat (v3.0.2).The corrected normalized data metric was applied to standard analysis.The top 50 variable genes were extracted for principal component analysis (PCA).The top 10 principal components were retained for UMAP visualization and clustering.Cell clustering was performed using the FindClusters function (resolution = 0.5) included in the Seurat R package.

Differential expression and survival prognosis analysis
The "Survival R package" was used to analyze differential expression, overall survival (OS), disease-specific survival, and progression-free survival of CD8 + T cell-associated genes based on single-cell sequencing and assays.Furthermore, we established the correlation between key genes and clinicopathological features of KIRC and constructed a prognostic nomogram and calibration curve in the TCGA-KIRC cohort.

Consensus clustering
A cluster analysis was performed using the "Consensus Cluster Plus R Package" [20] and agglomerative PAM clustering with a 1-Pearson correlation distances and resampling 80% of the samples for 10 repetitions.The optimal number of clusters was determined using the empirical cumulative distribution function plot.

Identification of differentially expressed genes (DEGs)
Differential mRNA expression between subtypes was studied using the Limma package (version: 3.40.2).The screening criteria were adjusted as p < 0.05 and |fold change| > 2.

Functional enrichment analysis and gene set enrichment analysis (GSEA)
The "clusterProfiler" package was used to evaluate the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of different CD8 + T cellassociated subtypes, with q-value and p-value thresholds of < 0.05.Moreover, the difference in gene sets between high and low expression of CD8 + T cell-associated genes was further assessed using the gene set enrichment analysis (GSEA) software (http:// www.broad insti tute.org/ gsea/ index.jsp).

Somatic mutation and immune landscape analysis
The KIRC somatic mutation data were obtained from the TCGA GDC data portal, and waterfall plots were created using the "Maftools" package.The relative percentages of different immune cell types were determined using the ESTIMATE and xCell algorithms, and the relative percentages of immune cell types between the two subtypes were compared by landscape plots.Moreover, TISIDB (http:// cis.hku.hk/ TISIDB/ index.php) was used to further evaluate the correlation between two subtypes and different immune indicators, including lymphocytes, immune inhibitors, immunostimulators, major histocompatibility (MHC) molecules, chemokines, and chemokine receptors.

Construction of CD8 + T cell-associated risk signature
Regression analysis was performed using the Lasso-Cox method via the glmnet package.A tenfold crossvalidation was set up to obtain the optimal model.The prognostic significance of genes involved was assessed using the Cox method, and the relationship between different risk scores and patient follow-up time, events, and changes in the expression of individual genes was analyzed.

Cell lines and cell culture
Human renal cancer cell lines (786-O and ACHN) were cultured in Dulbecco's modified Eagle medium (DMEM; Gibco Thermo Fisher Scientific, USA), containing 10% fetal bovine serum (Lonsera, Uruguay), and 1% penicillin-streptomycin solution (Keygen, China).All cell lines were purchased from the Shanghai Institutes for Biological Sciences and incubated in 95% humidified air at 37 °C and 5% CO 2 .

RNA extraction and RT-PCR
The RNA was extracted using the RNA extraction kit (Takara Kusatsu, Japan), and the Hiscript II First-Strand cDNA Synthesis Kit was used to synthesize the complementary DNA (Vazyme, China).Reverse transcriptasepolymerase chain reaction (RT-PCR) was performed using the MonAmpTM SYBR Green qPCR Mix (Monad Biotech, China).

Small interfering RNA
The small interfering RNAs (siRNAs) against the LAG3 gene were designed and synthesized by GenePharma Co. (China).

Cell proliferation and colony formation assays
For the cell proliferation assay, 1000 cells were seeded into 96-well plates for 0 h to 120 h, and 10 µL of the cell counting kit-8 (Keygen, China) solution was added per well.After a 2 h incubation at 37 °C, optical density at 450 nm (OD 450 nm) was measured using a microplate reader (Bio-Tek, USA).For the colony formation assay, cells were seeded into 6-well plates at a density of 1-2 × 10 [3] cells/well and incubated for 10 to 14 days at 37 °C.Next, the cells were washed using phosphatebuffered saline, fixed with 4% polyformaldehyde (Service Bio, China), and stained with 0.1% crystal violet solution (Keygen, China).Colonies containing > 50 cells were counted using the ImageJ 2X software 2.1.4.7 (Rawak Software Inc., Germany).

Wound healing and transwell assay
For the wound healing assay, cells were inoculated into 6-well plates and treated with si-/nc-LAG3.A straight scratch was created on the plate with a sterilized needle tip when the cell density was approximately 70%.The cell wound edge was marked and imaged under a microscope at the starting time point.After 0 h to 24 h, the migrated distance was measured and the wound closure percentage was calculated.For transwell assays, cells were inoculated into a 24-well transwell cell apical chamber containing the matrix gel (BD, USA) for evaluating invasion and gel-free for migration.The bottom and upper chambers contained the RPMI medium and fetal bovine serum-free medium, respectively.Cells invading the bottom chambers were fixed with 4% polyformaldehyde, stained with 0.1% crystal violet solution, counted, and imaged under a microscope.

Tissue samples and tissue microarrays
Formalin-fixed and paraffin-embedded prostate cancer tissue samples were collected from patients who underwent radical nephrectomy in the affiliated Zhongda Hospital of Southeast University, China, from April 2020 to November 2021.The study samples were from patients with KIRC, and the pathological diagnosis was confirmed by at least two pathologists.With the tumor as the center, normal tissues adjacent to the tumor were used as study materials, and two pairs of tissue microarrays were created with a 0.6 mm diameter.

Immunohistochemistry
Formalin-fixed and paraffin-embedded tissue was dewaxed and dehydrated using xylene and serially-diluted ethanol.The tissue sections were incubated at 121 ℃ in an autoclave for 5 min to extract the antigen, following which these were incubated with anti-LAG3-monoclonal antibody at 4 ℃ overnight, and the bound antibody (Proteintech) was incubated at 37 ℃ for 30 min.

Statistical analysis
The statistical analysis was performed using the R software (version 4.0.2).Multivariate Cox regression analyses were used to evaluate the prognostic significance.When p < 0.05 or log-rank p < 0.05, the difference was significant.

Identification of CD8 + T cell-associated gene subgroups by single-cell analysis
We used the GEO database to screen a single-cell sequencing dataset about KIRC GSE121636 [18], including subsets GSE3440844, GSE3440845, and GSE3440846, as well as the normal kidney single-cell sequencing dataset GSE131685 [18], including subsets, namely, GSE4145204 and GSEE4145205.The correlation between the pre-filtered (Fig. 1B) and post-filtered (Fig. 1C) base data was analyzed by excluding cells with less than 200 or more than 2500 genes detected and cells with mitochondrial content > 10% (Fig. 1A).The analysis suggested that the cells were more highly active and that the sequencing depth was not near saturation after data filtration.Gene clusters with a high degree of variation were obtained by the PCA analysis and displayed PCA plots of PC1 and PC2 after combining to remove the batch effects (Fig. 1D), as well as obtain the top 10 highly variable features (HVGs) in 2000 HVGs including S100A8, S100A9, JCHAIN, APOC1, TPSAB1, DEFB1, TPSB2, C1QB, C1QC, and SEPP1 (Fig. 1E).
The HVGs obtained above were included in the PCA analysis and 50 highly variable gene clusters were obtained using the ElbowPlot function and visualized by JackStrawPlot (Fig. 2A).Subsequently, KIRC cells and normal kidney tissue cells were divided into 14 clusters using the UMAP dimensionality degradation method (Fig. 2B).The top 10 genes in 14 clusters (Additional file 1: Table S1) were automatically matched to each cell type in the kidney tissue [21], ultimately obtaining eight different cell subtypes including the epithelium, CD8 + T cells, B cells, mast cells, NKT cells, neutrophils, Tregs, and plasma cells (Fig. 2C).A comparison of the original UMAP cluster with the artificially annotated cell cluster (Fig. 2B-C) revealed that B cells, NKT cells, neutrophils, and plasma cells were differentially expressed in KIRC and normal kidney tissue cells, with NKT cells being relatively highly expressed in KIRC (Fig. 2D).This analysis further revealed that the epithelial cells and mast cells were specific to the normal kidney tissue, whereas Tregs and CD8 + T cells were specific to KIRC, with CD8 + T cells exhibiting the highest expression (Fig. 2E).Moreover, different metabolic pathways between KIRC and normal kidney tissue were analyzed, with the KIRC tissue mostly enriched in DNA_replication, Tumor_prolif-eration_signature, Tumor_inflammation_signature, and IL-10_anti-inflammation_signaling_pathway, whereas the normal kidney tissue was enriched in MYC_targets, ECM-related_genes, and Cellular_response_to_hyposis (Fig. 2F).The artificially annotated cell clusters were used to study the metabolic pathways of each cell, with CD8 + T cell metabolic functions mostly enriched in Tumor_ proliferation_signature, Tumor_inflammation_signature, G2M_checkpoint, and DNA_replication (Fig. 2G), whose metabolic functions overlapped highly with those of KIRC.The above analysis indicated that CD8 + T cells could be an essential immune cell population in the formation, progression, and immune response to KIRC.

Identification of CD8 + T cell-associated subtypes by consensus clustering
The above single-cell sequencing analysis revealed that the top 10 CD8 + T cell-associated genes differentially expressed in KIRC were GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 (Additional file 1: Table S1).Subsequently, we verified the differential expression of these genes in the TCGA-KIRC cohort by heat map (Fig. 3A), paired differential expression analysis (Fig. 3B), and unpaired differential expression analysis (Fig. 3C).The results suggested that their expression was statistically significant in KIRC versus paracancerous tissues.Among these, GZMK, CD27, CCL4L2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 were highly expressed in KIRC and were considered oncogenes, whereas FXYD2 was lowly expressed.The differential expression of the above genes in cancer and normal tissues was further validated using the HPA database (https:// www.prote inatl as.org), which was consistent with the TCGA-KIRC cohort analysis.However, the database did not have immunohistochemical plots for CCL4L2, DUSP4, and TRBV20-1 (Additional file 3: Figure S1, Additional file 2: Table S2).Moreover, the number of clusters with the highest average consistency within the group was K = 2 (Fig. 3D), and the distribution curve was the greatest at K = 2 (Fig. 3E).Therefore, K = 2 was selected to perform clustering, and the consistent cluster (K = 2) is depicted in Fig. 3F, with 231 KIRC patients in subtype C1 and 276 KIRC patients in subtype C2.Differential expression and OS analyses of the top 10 CD8 + T cell-associated genes in C1 and C2 subtypes further revealed that GZMK, CD27, CCL4L2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1 were highly expressed in the C2 subtype, whereas FXYD2 was highly expressed in the C1 subtype (Fig. 3G).The OS was shorter in the C2 subtype compared to the C1 subtype (Fig. 3H).These results were similar to the differential expression analysis of the TCGA-KIRC cohort, indicating that the C2 subtype was the oncogenic group for KIRC.

Identification of differentially expressed genes and signal pathways in different CD8 + T cell-associated subtypes
The above analysis revealed a poor prognosis in the oncogenic subtype (C2), whereas the low expression subtype (C1) of CD8 + T cell-associated genes had a better prognosis.Consequently, we identified key DEGs and signaling pathways in each subtype to understand the molecular mechanisms.A total of 507 abnormally regulated genes were identified, including 15 up-regulated and 492 down-regulated genes (Fig. 4A-B).Because of a few up-regulated genes, we only performed the KEGG analysis on down-regulated genes and identified them to be mostly enriched in cytokine-cytokine receptor interaction, primary immunodeficiency, T cell receptor signaling pathway, cell adhesion molecules (CAMs), antigen processing and presentation, and chemokine signaling pathway (Fig. 4C).The up-regulated genes were largely enriched in glutamate and leukotriene activities, including hypoglycin A gamma-glutamyl transpeptidase activity, leukotriene C4 gamma-glutamyl transferase activity, and glutathione hydrolase activity (Fig. 4D).However, the down-regulated genes majorly corresponded to immune-related activities, including immune system processes, immune responses, and adaptive immune responses (Fig. 4E).Moreover, the GSEA analysis further suggested that C1 and C2 subtypes significantly differed in gene set enrichment in immune cells (T cell receptor signaling pathway, antigen processing and presentation, natural killer cell-mediated cytotoxicity, B cell receptor signaling pathway, and leukocyte transendothelial migration) (Fig. 4F) and oncogenic pathways (VEGF signaling pathway, P53 signaling pathway, DNA replication, apoptosis, and pathways in cancer) (Fig. 4G).These results suggested significant differences in gene expression and related pathway enrichment between C1 and C2 subtypes, which could be the underlying mechanism leading to the different prognosis of KIRC.
Because the study was based on CD8 + T cells, which could potentially influence the activation or silencing of tumor immune responses, we analyzed the tumor immune microenvironment in two subtypes.The StromaScore, ImmuneScore, ESTIMATEScore, and Microenviron-mentScore were higher in the C2 subtype compared to the C1 subtype (Fig. 5A-C, 5F-G) although the xCell algorithm's StromaScore subgroup type was not meaningful (Fig. 5E).Next, the infiltration differential map of 64 immune cell types was assessed by xCell with the LM64 signature matrix (Fig. 5D).The quantitative analysis indicated that the aDC, CD8 + T cells, CD8 + Tcm, macrophages, NKTs, Th1 cells, and Th2 cells in the C2 subtype were significantly higher than those in the C1 subtype, whereas monocytes and smooth muscle cells were lower than those in the C1 subtype (Fig. 5H).Subsequently, a comprehensive analysis of the immune gene correlation between different subtypes revealed that the expression of both immune checkpoint-related genes and MHC molecules was significantly higher in the C2 subtype than in the C1 subtype (Fig. 5I-J).For immunostimulators and chemokines, the expression of CD27, CD48, ICOS, CCL5, CXCL9, and CXCL10 in the C2 subtype was significantly higher than in the C1 subtype (Fig. 5K-L).

Construction and validation of CD8 + T cell-associated risk signature
We used the LASSO regression analysis and Cox univariate analysis to construct a risk prognostic model based on CD8 + T cell-associated genes.After integrating  6A-B).Subsequently, it was finally determined by the Cox univariate analysis that GZMK, FXYD2, LAG3, RGS1, and DUSP4 displayed prognostic significance (Fig. 6C).Next, the relationship between survival status and the five genes above was investigated.Our results demonstrated that the number of survival states was considerably higher in the anti-cancer cohort (C1 subtype) compared to the oncogenic cohort (C2 subtype) (Fig. 6D).The risk scores for the above five genes were further quantified using the KM analysis in the TCGA-KIRC cohort, where high-risk scores corresponded to poorer OS (Fig. 6E) and 1, 3, and 5-year survival rates of 0.67, 0.67, and 0.69, respectively (Fig. 6F).

Identification of LAG3 as the most critical CD8 + T cell-associated gene in KIRC
The correlation of CD8 + T cell-related genes (GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1) was further explored with prognosis and diagnosis in KIRC.We further validated the Fig. 5 Immune landscape of C1 and C2 subtypes of CD8 + T cell-associated genes.Subtypes were assessed using the ESTIMATE algorithm for stromalScore (A), immueScore (B), and ESTIMATE score (C).Immune cell infiltration (D) was assessed using the xCell algorithm, as well as stromalScore (E), immueScore (F), and microenvironment Score (G) for different subtypes.Differential expression of genes related to immune cells (H), immunosuppression (I), MHC molecules (J), immune enhancement (K), and chemokines (L) between C1 and C2 subtypes.C1 was a low-expression subtype of CD8 + T cell-associated genes and C2 was a high-expression subtype OS, disease-specific survival, progression-free interval, diagnostic receiver operating characteristic (ROC), and time-dependent ROC curves of the above 10 genes in the TCGA-KIRC cohort.The results suggested that the area under the ROC curve was > 0.85 for all genes, except FXYD2, and the diagnostic effect was good (Additional file 4: Figure S2A-J).LAG3 exerted a significant pro-carcinogenic effect and KIRC patients with high LAG3 expression had shorter OS (p = 0.008), diseasespecific survival (p = 0.006), and progression-free interval (p = 0.043) (Additional file 4: Figure S2H).In contrast, FXYD2 and DUSP4 displayed a significant anti-carcinogenic effect with all p < 0.05 (Additional file 4: Figure S2E-F), whereas the other single genes exerted no great prognostic significance (Additional file 4: Figure S2).Notably, only LAG3 exhibited a relatively good timedependent curve, AUC 1-year = 0.592, AUC 3-year = 0.574, and AUC 5-year = 0.580 (Additional file 4: Figure S2H).Further analysis of its association with CD8 + T cells revealed that the expression of LAG3 and DUSP4 was positively correlated with CD8 + T cell enrichment (Fig. 7A-B), whereas FXYD2 was not (Fig. 7C).The correlation between key genes and immunotherapy was further explored by TMB scores, suggesting that LAG3 was positively correlated with TMB scores (ρ Spearman = 0.14 and p = 0.006), whereas DUSP4 was insignificant (Fig. 7D,  F).The study revealed that the expression of LAG3 and DUSP4 was not associated with TIDE scores (Fig. 7E, G).Compared to the expression of DUSP4, LAG3 displayed a higher correlation with CD8 + T cells and a function in immunotherapy.Therefore, we followed up with an indepth exploration of LAG3.We further verified the correlation between LAG3 expression and different T-cell subtypes in the TCGA-KIRC disease cohort.We found that the LAG3 level was positively correlated with several T cells, including CD8 T cells, cytotoxic cells, T cells, T helper cells, Tcm, Tem, TFH, Th1 cells, Th2 cells, and Tregs, and negatively correlated with Th17 cells, with no correlation with Tgd (Additional file 5: Figure S3).Therefore, LAG3 was identified as the most critical T Fig. 6 Construction and validation of CD8 + T cell-associated genetic risk prognostic model in the TCGA-KIRC cohort.A-B Correlation of 10 CD8 + T cell-associated genes with OS using the Lasso-Cox analysis.C The prognostic value of genes for OS was assessed using the univariate Cox analysis.D Risk score distribution, survival status, and heat map of five valuable genes.Kaplan-Meier plots (E) and ROC plots (F) of CD8 + T cell-associated genetic risk prognostic models cell-related prognostic gene in KIRC.Further investigation of the association of LAG3 expression with sunitinib and sorafenib sensitivity indicated that the high LAG3 expression group (G2) had lower sunitinib IC 50 and higher sorafenib IC 50 (Fig. 7H-I).Moreover, LAG3 was correlated with multiple immune-related genes (Fig. 7K).This study revealed that the LAG3 expression was not associated with stemness scores between the high and low LAG3 expression groups (Fig. 7J).
The single-cell analysis revealed that LAG3 expression was mostly concentrated in zone 1 with scattered distribution in zones 2, 5, 6, 7, and 10 (Fig. 8A).Cluster identification demonstrated that LAG3 was concentrated in CD8 + T cells in KIRC, whereas it was highly scattered in NKT or B cells in both normal kidney tissues and KIRC (Fig. 8B).The GSEA pathway after single-gene differential analysis suggested that its pathway included the PID_CD8_TCR_PATHWAY, WP_INTER-ACTIONS_BET WEEN_IMMUNE_CELLS_AND_ MICRORNAS_IN_TUMOR_MICROENVIRONMENT, W P _ C A N C E R _ I M M U N OT H E R A P Y _ BY _ P D 1 _ BLOCKADE and REACTOME_IMMUNOREGULA-TORY_INTERACTIONS_BETWEEN_A_LYMPHOID_ AND_A_NON_LYMPHOID_CELL, whose NES > 1.9, p.adj < 0.05, and FDR < 0.05 (Fig. 8C-F).In addition, the combined enrichment scores of each sample on multiple pathways revealed the most significant positive correlation between LAG3 and Tumor_Inflammation_Signature, apoptosis, and IL-10_Anti-inflammatory_Signaling_ Pathway (Fig. 8G-I).

Validation of LAG3 as a key CD8 + T cell-associated gene in KIRC in in vitro experiments
An analysis of the CCLE dataset revealed a significantly high expression of LAG3 in A-704, OS-RC-2, SNU-349, KMRC-3, SLR 26, and ACHN kidney cancer cell lines (Additional file 6: Figure S4).Next, the RT-PCR results of existing kidney cancer cell lines in the laboratory found LAG3 to be highly expressed in 786-O and ACHN cell lines, which were selected for subsequent validation (Fig. 9A).The knockdown efficiency of LAG3 small interfering reagents in 786-O and ACHN cell lines by RT-PCR was validated (Fig. 9B).Subsequently, the growth curve of the siLAG3 group was significantly slower in the 786-O and ACHN cell lines than that of the control group in the CCK8 experiment (Fig. 9F).Compared with the control group, the number of cell clones in the siLAG3 group was significantly lower in 786-O and ACHN cell lines by clone experiments (Fig. 9C).The scratch experiment revealed that the siLAG3 group was less efficient at healing scratches at 24 h than the control group (Fig. 9D).Cell migration and invasion were detected to be significantly down-regulated in the siLAG3 group by transwell assay (Fig. 9E).Combined with the preliminary results of LAG3 differential expression in KIRC versus normal kidney tissues in the HPA database (Additional file 3: Figure

Discussion
Kidney cancer constitutes the 14th most common malignancy worldwide with 431,288 new diagnosed cases and 179,368 new deaths in 2020 [3].The major etiologies of kidney cancer include hypertension, obesity, and smoking [22].The ICI-based combination therapy demonstrated excellent clinical efficacy in several large clinical trials and is now the first-line care standard for patients with advanced or metastatic renal cancer with a low OS at first diagnosis [23][24][25][26].Despite the pivotal function of the PD-1/CTLA-4 axis in the treatment of RCC has greatly improved clinical outcomes compared to previous treatment options, the majority of patients with RCC did not achieve durable clinical benefits after ICI-based combination therapy [23][24][25][26].Therefore, it is highly essential to investigate the tumor immune microenvironment and explore novel immunotherapeutic targets, and ultimately optimize systemic treatment for RCC.
In response to the immune microenvironment of renal cancer, Finke [6] and Stein [15] provided reviews on the immunology and immunobiology of renal cancer, respectively.In most solid tumors, the degree of CD8 + T cell infiltration was positively correlated with good prognosis for tumor patients [27].CD8 + T cells exerted a direct cytotoxic effect on target cells and performed a critical role in anti-cancer immunity.However, the expression of suppressive molecules (PD-1 and CTLA-4) on the surface of CD8 + T cells increases in response to sustained stimulation by tumor-specific antigens, and their function decreases and eventually reaches the exhausted state, as demonstrated in multiple cancer models [13,14,28].Therefore, blocking their inhibitory signaling using anti-PD1/CTLA-4 antibodies could rejuvenate exhausted CD8 + T cells, enhance their cytotoxicity, promote tumor cell lysis, and restrict tumor metastasis [29].However, kidney cancer had a distinctive immune profile such that the degree of CD8 + T cell infiltration is positively correlated with poor prognosis, and the specific mechanism was unclear.Several hypotheses have been proposed to explain this paradoxical phenomenon.First, the activation status and virulence potential of CD8 + T cells were highly specific in kidney cancer, where stem-like TCF1 + or PD-1 + TIM3 − LAG3 − CD8 + T cell subsets contribute to the anti-cancer immune effect [30][31][32].Second, the low density of tertiary lymphoid structures generated numerous immature DC cells, causing the infiltration of polyclonal CD8 + T cells that could not recognize tumor-associated antigens [11,27,32].Third, the absence of tumor-specific genes, such as the relative absence of PBRM1 mutations in highly CD8 + T cells RCC, which was often associated with a better prognosis [33].Finally, specific metabolic dysregulation of CD8 + T cells in RCC restricted CD8 + T cell activation and did not recover through the PD-1 axis inhibition [34].
The high significance of CD8 + T cells in the immunotherapy of kidney cancer and the results of single-cell sequencing analysis (CD8 + T cells were significantly differentially expressed in KIRC and kidney tissue) can be used to construct a prognostic model of CD8 + T cellassociated genes to guide clinical decisions in KIRC.First, the top 10 DEGs in CD8 + T cells of KIRC were obtained by cluster identification in single-cell analysis, including GZMK, CD27, CCL4L2, FXYD2, LAG3, RGS1, CST7, DUSP4, CD8A, and TRBV20-1.Subsequently, the above gene expression was divided into two subtypes by cluster analysis, that is, CD8 + T cell-associated gene low expression (C1) and high expression (C2) subtypes.The cluster analysis results were used for grouping, and the DEGs, pathway enrichment, and mutated genes between C1 and C2 subtypes were comprehensively studied.The up-regulated genes were largely enriched for glutamate and leukotriene activity, whereas the down-regulated genes were enriched for immune-related activities.Furthermore, VHL and STED2, the most commonly mutated genes in primary KIRC, displayed remarkably high mutation rates in the C2 subtype compared to the C1 subtype, predicting a poor prognosis such as metastasis or drug resistance in the C2 subtype [35,36].We assessed the KIRC immune microenvironment in different subtypes and revealed that the stromal and immune scores were significantly higher in the CD8 + T cell-associated gene high expression subtype than in the C1 subtype, with significant differences between the two subtypes in immune infiltrating cells and immune-related molecules.Therefore, LASSO regression and Cox univariate analyses were used to construct the CD8 + T cell-associated risk prognostic model: RiskScore = − 0.291858656434841*GZM K − 0.192758342489394*FXYD2 + 0.625023643446193* LAG3 + 0.161324477181591*RGS1 − 0.38016904532889 5*DUSP4 − 0.107221347575037*TRBV20-1, which will assist the clinicians to assess the prognosis for survival, immune status, and drug selection in patients with KIRC.Relevant studies were investigated to identify CD8 + T cell-related genes in kidney cancer.Genomics, radiology, and artificial intelligence modalities can be used to identify renal cancer differentiation more easily and earlier, and predict its Fuhrman grade and responsiveness to immunotherapy response, thus assisting the clinicians in defining the risk of stratification of the disease, treatment choices, follow-up strategies, and prognosis [42,43].
Using tumor single-gene related studies as the reference [44], we performed single-gene analysis of the above genes and finally identified LAG3 as the most valuable CD8 + T cell-associated gene in KIRC.Lymphocyte Fig. 10 LAG3 is highly expressed in the KIRC tissue activating gene 3 (LAG3) or CD223 is highly expressed in different T cells, CD8 + T cells, CD4 + T cells, and Tregs, to maintain homeostasis [37].However, persistent tumor-associated antigen stimulation causes its chronic expression, ultimately promoting T-cell exhaustion in cancers [37,38].Therefore, LAG3 serves as the third clinical checkpoint in case of limited or even no response in 60 to 80% of cancer patients in PD1/CTLA-4 axis immunotherapy [39].Currently, several clinical trials are being conducted on immunotherapies targeting LAG3 in combination with PD1/CTLA-4 axis inhibitors to treat cancers including kidney cancer 40, 41.In the current study, we analyzed for the first time the distribution of LAG3 in renal cancer and renal tissue using single-cell analysis and investigated its expression, prognostic significance, immune microenvironmental relevance, and pathway enrichment using bioinformatics techniques.We confirmed that LAG3 promotes the progression and metastasis of renal cell carcinoma and is positively correlated with CD8 + T cells using cell phenotype studies and immunohistochemistry. Presently, immunotherapy targeting LAG3 is largely used for melanoma, pancreatic cancer, and hematological tumors, with only a few studies on renal cancer.We elucidated the function of LAG3 in KIRC.We believe our findings will provide a preliminary basis and direction for LAG3-targeted immunotherapy and even CAR-T therapy in patients with kidney cancer.
The present study had certain limitations.First, heterogeneity obtained in retrospective studies needs to be verified by conducting prospective studies.Second, we only applied the top 10 CD8 + T cell-associated genes to construct the risk prognosis model, which lacked comprehensiveness and extensiveness.Third, we only validated the function of LAG3 in KIRC at the in vitro level and lacked in vivo experiments, as well as the expression mechanism and CD8 + T cell activity warrant further exploration.In conclusion, more basic and large clinical trials are required to validate these findings (Table 1).

Fig. 1
Fig. 1 Pre-procession of single-cell sequencing data.A Genetic data, molecule numbers, and mitochondrial genome percentage in each cell.B Correlation between the primary data before filtration.C Correlation between primary data after filtration.D RPCA plot after the removal of batch effects.E Highly variable genes

Fig. 2
Fig. 2 Identification of sub-groups in single-cell analysis in KIRC.A Determination of the principal components using ElbowPlot and JackStrawPlot() functions.B Cell cluster plotting using UMAP.C Cell cluster plotting after manual annotation of the top 10 marker genes.Plots of the distribution of shared cells (D) and unique cells (E) in normal kidney tissues and kidney cancer.F Metabolic pathways of normal kidney tissue and kidney cancer.G Metabolic pathways of different immune cells

Fig. 3
Fig. 3 Identification of CD8 + T cell-associated subtypes by consensus clustering.A Gene expression profiles of 10 CD8 + T cell-associated genes in the TCGA-KIRC cohort.Paired (B) and unpaired (C) expression analyses of 10 CD8 + T cell-associated genes.Sample cluster consistency (D) and area under the distribution curve (E) for K from 2 to 10. F Consistent cluster (K = 2) of 10 genes.G Expression profiles of 10 CD8 + T cell-associated genes in different subtypes.H Kaplan-Meier curves for OS in C1 and C2 subtypes.C1 was a low-expression subtype of CD8 + T cell-associated genes and C2 was a high-expression subtype

Fig. 4
Fig. 4 Identification of differentially expressed genes, their potential signaling pathways, and somatic mutations in C1 and C2 subtypes.DEG expression profiles (A) and volcano plots (B) of different subtypes.C KEGG analysis of down-regulated DEGs.GO analysis of up-regulated (D) and down-regulated (E) DEGs.Potential signaling pathways between C1 (F) and C2 (G) subtypes.Visualization of the top 10 mutated genes in the C1 (H) and C2 (I) subtypes of the CD8 + T cell-associated genes.C1 was a low-expression subtype of CD8 + T cell-associated genes and C2 was a high-expression subtype

Fig. 7
Fig. 7 Single-gene analysis of LAG3 in CD8 + T cell-associated genes.The expression and correlation with CD8 + T cells for LAG3 (A), DUSP4 (B), and FXYD2 (C) in the TCGA-KIRC cohort.Correlation of LAG3 with tumor mutation burden (D) and responsiveness to immune checkpoint inhibitors (E).Correlation of DUSP4 with tumor mutation burden (F) and responsiveness to immune checkpoint inhibitors (G).IC 50 of sunitinib (H) and sorafenib (I) in different LAG3 expression subgroups.(J) LAG3 stemness score.(K) Relevance of LAG3 to immune-related genes

Fig. 8
Fig. 8 Cell population expression and pathway analysis of LAG3 in single-cell analysis.A LAG3 UMAP plot in CD8 + T cell-associated subgroups.B Expression of LAG3 in different cell populations in normal kidney and renal cancer tissues.C-F Analysis of potential signaling pathway of LAG3 by GESA.G-I Analysis of the enrichment of LAG3 using the GSVA algorithm

Fig. 9
Fig. 9 Validation of LAG3 as a key oncogene in KIRC in in vitro experiments.A Validation of LAG3 expression in renal cancer cell lines.B Knock-down efficiency of LAG3, respectively, in 786-O and ACHN cell lines.C LAG3 clone tests in 786-O and ACHN cell lines.D LAG3 scratch tests in 786-O and ACHN cell lines.E LAG3 migration and invasion assays in 786-O and ACHN cell lines using a 24-well plate.F LAG3 CCK8 assay in 786-O and ACHN cell lines